Moving boundary approximation for curved streamer ionization fronts: Numerical 

tests 



On 
O 
O 
(N 

G 

03 



43 

a: 



c/3 

11: 

o ■ 

>V 
43 . 
Oh. 



> 
On 



O 

On 
O 



X 



Fabian Brau 1 , Alejandro Luque 1 , Benny Davidovitch 2 and Ute Ebert 1,3 
1 Centrum Wiskunde & Informatica (CWI), P.O.Box 94079, 1090 GB Amsterdam, The Netherlands 
2 Physics Department, University of Massachusetts, Amherst MA 01003 and 
3 Dept. Physics, Eindhoven Univ. Techn., The Netherlands 
(Dated: January 14, 2009) 

Recently a moving boundary approximation for the minimal model for negative streamer ioniza- 
tion fronts was extended with effects due to front curvature; this was done through a systematic 
solvability analysis. A central prediction of this analysis is the existence of a nonvanishing electric 
field in the streamer interior, whose value is proportional to the front curvature. In this paper we 
compare this result and other predictions of the solvability analysis with numerical simulations of 
the minimal model. 



I. INTRODUCTION 

Streamers characterize the initial stages of electric 
breakdown in sparks, lightning and sprite discharges; 
they occur equally in technical and natural processes 
0, 0, HI- They are growing plasma channels that ap- 
pear when strong electric fields are applied to ioniz- 
able matter. The essential features of negative (anode- 
directed) streamers in a non-attaching gas such as argon 
or nitrogen can be described by the so-called minimal 
model H, |, 1, 0, I, & M EH, Gl El G3- This model 
consists of a set of three coupled partial differential equa- 
tions for the electron density a, the ion density p and the 
electric field E. In dimensionless units the model reads 

d t a- V- (crE) - DV 2 cr = cr|E| a(|E|), (1) 
d tP = a\E\ a(|E|), (2) 
V-E = p-cr, E = -V>, (3) 

where D is the electron diffusion coefficient and where 



|E|) = 



- P -V|E| 



(4) 



A general discussion of the physical dimensions for this 
model can be found, e.g., in @, 0, 0, EH- The model 
is based on a continuum approximation with local field- 
dependent impact ionization reaction. Equations ([T]) and 
@ are the continuity equations for the electrons and the 
ions, taken as immobile due to their much larger mass, 
while Eq. is the Coulomb equation for the electric 
field generated by the space charge p — a of electrons 
and ions. Although discharges in air require extensions 
of the model, simulation results of negative air stream- 
ers f req uently resemble the minimal model remarkably 

well diEi]. 

Many simulations @, 0, 0, 0, EI EH, El, El have 
shown that streamers form a thin curved space charge 
layer which separates the ionized interior region, Q~ , 
from the nonionized exterior region, il + . This narrow 
charged layer (the ionization front) enhances the electric 
field in £l + ahead of the front and screens it partially 
in f2~ . In strong background fields after some transient 
evolution, the width of the ionization front can be much 
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FIG. 1: (Color online) On the left: representative solution 
(net charge density) of the minimal PDE model with cur- 
vature K and width (see Eq. ([9])) of the ionization front. 
The electric field E is pointing downward and the negative 
front propagates upward. On the right: depiction of the cor- 
responding moving boundary approximation (MBA) with the 
ionized region, the non- ionized region, Q + , and the sharp 
interface, F. 



smaller than its radius of cur vature 0, E3, El- This 
separation of scales enables one to consider the front as 
an infinitesimally thin, sharp moving interface T(t). In 
Fig. [TJ wc show a representative snapshot of net charge 
density of the minimal model ((TJ) — (|SJ) that shows the sep- 
aration of scales, and depict the corresponding moving 
boundary approximation. The original nonlinear dynam- 
ics is then replaced by a set of linear field equations (fre- 
quently of diffusive or Laplace type) on both sides of T(t); 
the regions on both sides of T(t) are denoted as tt + and 
f2~. The linear fields in these regions are determined 
by boundary conditions on both sides of the interface, 
r(£) + , r(t)~, respectively, and on the outer boundaries 
(assumed to be located far away from T(t)); the nonlin- 
earity enters through the motion of the boundary. The 
interface dynamics is typically related to gradients of the 
Laplacian fields in its vicinity. 
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In the context of streamer dynamics, the concept of an 
interfacial approximation was probably first sketched by 
Sammer in 1933 [17j |; later it was developed further by 
Lozansky and Firsov in the Russian literature and in En- 
glish in [18J. They considered the streamer interior, f2 _ , 
as ideally conducting, i.e., the electric potential <p to be 
constant in the interior. The exterior, is nonionizcd 
and therefore does not contain space charges; the electric 
potential here solves 



VU = 



fi+. 



(5) 



The interface was assumed to move with the local elec- 
tron drift velocity 



v = V</> 4 



(0) 



Hence, superscripts ± attached to fields, potential and 
densities indicate their limit values as they approach the 
interface from f2 + and respectively. In particular, 
we denote <j) + = 4>\y+ and <p~ = (f>\r- ■ 

However, this simplest moving boundary approxima- 
tion is mathematically ill posed; in the context of similar 
models in fluid dynamics, this is explained for example in 
Ref. [IH and references therein. To resolve this problem, 
the boundary condition </> + — (f>~ = was replaced by the 
regularizing boundary condition 



Qo(n • 



|n-V0 + |>l 



V(/> H 



(7) 



which was proposed in [201 ] and derived in planar front 
approximation in [2~D ] . The boundary condition accounts 
for the finite width of the charged layer that leads to a 
finite variation of the electric potential across the front. 
The boundary condition in the limit of large electric fields 
actually turns out to be identical to the "kinetic under- 
cooling" boundary condition that was applied to crystal 
growth under certain conditions 1221. 231 . Solutions of 
the model ©-CO) are discussed in [20, [24|, and the anal- 
ysis in [25[ shows that the boundary condition ((7|) indeed 
regularizes the problem. This moving boundary approxi- 
mation is com par ed with solutions of the minimal model 

©-© in MlM- 

In a recent paper (27[, effects associated with curva- 
ture of the front were considered. The moving bound- 
ary conditions for a slightly curved front dynamics were 
systematically derived from the original nonlinear field 
equations |T|)-([3]), with D = 0, using the following proce- 
dure. A perturbation of a planar front is assumed whose 
curvature in the direction transverse to the front motion 
is much smaller than the front width, 



«1, 



(8) 



where £~ and k are respectively the width and the cur- 
vature of the front respectively. The width of the front, 
l~ , is taken as the decay length in the ionized region of 
the net charge density of the planar front with D = 
and reads 123 



E+ 



(9) 



where <r~ is the value of the electron density far behind 
the planar front whose expression is given by Eq. (|19|) . 
see for example Refs. [H, [13, HE]- £~(E + ) is a monoton- 
ically decreasing function of E + and tends to 1 (in our 
dimensionless units) for E + — > +oo. The computation 
is carried out in first order in e. Solvability analysis is 
then used to connect the perturbed values of the fields 
ahead and behind the curved front to derive the moving 
boundary approximation. Similarly to other well-studied 
problems such as solidification dynamics p9l . l30l | , this ex- 
pansion around the planar front solution is asymptotic 
and does not necessarily converge. However, such solv- 
ability analysis provides valuable approximation for the 
nonlinear dynamics of the propagating front as long as 
e remains sufficiently small. Furthermore, notice that 
such an analysis could not be performed on the streamer 
model with D > as the fronts are pulled [U HH, [33| . 
However, the leading edge that pulls the front is diffusive, 
and it is a physically and mathematically meaningful ap- 
proximation to neglect electron diffusion in strong fields, 
where electron motion due to drift dominates over the 
diffusive motion [28|, [34[ ■ 

The complete model derived in Ref. (2?| reads 



V z cj) = 
V 2 6 = 



in W 
in f2~ 



(10) 
(11) 



with the moving boundary conditions 



h-V<f>- = Q 2 (n-Vc/>+)K (12) 
+ -<T = Q (n- V0 + ) + Qi(n- V</> + )k (13) 
v n = n-V(f>+. (14) 

The coefficients Qi depend on the electrostatic field 
ahead of the front, and are given by analytic expressions 
derived from the planar front solution as follows 



/ y — % 
My) = / dz —, — \> 
Jo p(z,y) 

fy (y - x)a(x) f x (y-z) 

My) = -y / dx — — / dz— — - 

Jo xp[x,y) J p(z,y) 2 
1 f y , (y 2 - x 2 ) 



(15) 



dx- 



o (y) Jo ~~ xp(x, y) 2 



V 



dx 



y-x 



[p{x,y)y-a {y-x)] 

3 



y 



° (y) Jo p(x,y) (a- (y)) 2 ' 



Qa(tf) 



where 



° (y) : 



P(x,y) = / dpa{p), 

J\x\ 

ry 

<r~(y) = / dpa(n), 



(16) 
(17) 

(18) 
(19) 



with a(x) = e x / x . The quantity p is related to the ion 
density profile of the uniformly translating planar front 
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solution of the minimal model H])-© with D = 0, see 
for example Refs. @, H3, HI] . 

The boundary condition (fT2|) implies that the electric 
field just behind the ionization front is not completely 
screened but that it is proportional to the curvature. 
This implies that the ideal conductivity approximation 
in the streamer interior (<f> = in Q~) must be relaxed. 
In Ref . [27| , the streamer interior was therefore approxi- 
mated by assuming charge neutrality (V 2 <j> = in 
Consequently, the boundary condition (|12p introduces 
new physics and leads to a new type of moving boundary 
problem. 

The purpose of this work is to study the validity of this 
moving boundary model and, for the reason just men- 
tioned, especially the validity of the boundary condition 
(|12[) by comparing it to the results obtained from numer- 
ical simulations using the minimal model ©-(HI). No- 
tice that the numerical simulations are performed with 
a non- vanishing electron diffusion coefficient D. There 
are essentially two reasons for that. First, a vanishing 
diffusion coefficient leads to a continuum model with an 
electron discontinuity that cannot be simulated with the 
numerical methods developed for nonvanishing diffusion. 
Second, we want to test our boundary conditions on a 
realistic model and see if our moving boundary model is 
robust against some changes in the underlying minimal 
model. 

Another relation derived in Ref. (2?| , that does not ap- 
pear explicitly in the model pH)) - (rTi|) . can also be tested 
against numerical simulations. This is the curvature cor- 
rection to the value of the electron density behind the 
front 



where 



CTback = CT (n • V0 + ) + Q 3 (n • V0 H 



Qs(y) = y[ V dJ y - x)a(x) 
Jo 



(20) 



(21) 



xp(x,y) 

The paper is organized as follows. In Sec. HU we 
describe the method used for comparing the moving 
boundary approximation with the simulation data, and 
in Sec. IIIIl we describe in detail the results of our com- 
parison concerning the boundary conditions (TT2")) and (fT3"|) 
and also Eq. l[2"0"]). 



METHOD FOR COMPARING THE MOVING 
BOUNDARY APPROXIMATION WITH 
SIMULATIONS OF THE MINIMAL MODEL 



II. 



In this section, we test the boundary conditions (jT2|) 
and (Tl3|) as well as Eq. ([20]) against results of simulations 
of the minimal model (H}-© in two dimensions. The elec- 
tric field and the electron density behind the front are es- 
sentially constant over a significant interval, therefore it 
is relatively easy to extract their values from simulation 
data without introducing significant errors. The compar- 
ison with predictions of Eqs. (TT2")) and (|2T)|) allow to test 



the model with confidence. In contrast, as explained be- 
low, due to some arbitrariness of the precise location of 
4> + in the simulations and since the potential varies sig- 
nificantly over short distances, the comparison between 
Eq. (|13[) and the simulations is not quite conclusive. 

In order to test our boundary conditions, we need to 
evaluate the profiles of the net charge density, of electric 
field and potential and of the electron density along some 
given axis of the two-dimensional simulated streamer. In 
this paper we consider a streamer that evolves from ini- 
tial conditions with a mirror symmetry y —> —y. Fur- 
thermore, we restrict our analysis to the symmetry axis 
of the streamer, y = 0, since this allows easier evaluation 
of the curvature and the various fields. 



A. Numerical simulations 

The minimal PDE model ©-© (with D = 0.1 @, 
0, IE]) is solved numerically in two dimensions on adap- 
tively refined comoving grids with a second-order explicit 
Runge-Kutta time inte gra tion. The algorithm is de- 
scribed in detail in Ref. [1J| for three-dimensional, cylin- 
drically symmetrical geometries. It is trivially adapted to 
planar two-dimensional systems, as previously discussed 
in Refs. [2l], [26|. The highest spatial resolution in the 
area around the streamer head was Aa; = Ay =1/4 for 
all simulations. The simulation domain was < x < 2048 
and —1024 < y < 1024. The initial conditions were 
an electrically neutral Gaussian seed of width 16 and 
height 2.4 • 10~ 5 . We used four different values for the 
background electric field applied between the electrodes, 
namely E x = —0.5, —1, —1.5 and —2. The simulations 
are the same as in [2l| ; and for the actual density and 
field configurations, we refer to the figures in that paper. 



B. Extracting relevant quantities from the 
simulation data 



For each value of E 00l we collected at constant time 
steps, up to the time of branching, the values of the cur- 
vature of the front, of the enhanced electric field (defined 
as the maximum of the electric field along the symmetry 
axis) and the profiles of the electric potential and elec- 
tron density. These are the ingredients of the boundary 
conditions (12]) and (US]) and of Eq. (J20]) that we test in 
this paper. 

To test both relations (|12p and (Ti~3")) using a unique 
procedure, we consider Eq. (128) of Ref. [27[ (where the 
leading contribution of the planar front is added here), 

0(0) - 0(x) = Q {E + ) + Qi{E+)k - Q 2 {E+)kx, (22) 

where |x| 3> and where x = corresponds to the 
position of the tip of the front, i.e. the position of the 
discontinuity line T(t). This equation predicts that the 
correction to the potential profile due to curvature is a 
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FIG. 2: (Color online) Profiles of the electric potential and 
of the negative net charge density along the symmetry axis, 
y = 0, for Eoo — — 0.5 and t = 490. The linear regression for 
the linear part of the potential is also plotted. The slope of 
this linear part corresponds to the electric behind the front, 
£.-(sim) rp^g (jiff erence between the simulated value of the 
potential at x = and the value of the extrapolation of the 
linear part at the same location gives the jump in the electric 
potential. 

linear function of the variable x behind the front in an in- 
termediate region between the inner region and the outer 
region [27\. The position x = where the potential (f)(0) 
is evaluated, is taken as the location of the maximum of 
the negative net charge density, while E + is identified 
with the maximum of the electric field along the symme- 
try axis (since we restrict our analysis here to this axis). 
The slope of the linear curve (|2"2")l yields the boundary 
condition (fl"2"l) and expresses the electric field behind the 
front. The comparison between the slope in the simu- 
lations and the corresponding expression in Eq. (|22[) is 
then a direct test of Eq. (fT2"|) . Equation (f2"2"|) can also be 
used to test the relation (fT3")) . For this purpose, we ex- 
trapolate the linear part of the potential, 0n n , up to the 
tip of the front (x = 0) and compare the difference be- 
tween this value obtained for the potential, 0n n (O), and 
the value of the potential obtained from the simulation 
at x = 0, </>(0). Indeed we have 

0(0) - 0iin(O) = Qo(E+) + Qi(E + ) K , (23) 

where 0(0) can be measured from the simulated potential 
and 0ii n (O) is obtained from the linear regression; the 
procedure is illustrated for — —0.5 and t = 490 in 
Fig. [2j We then can compare the simulated potential 
jump with the theoretical value on the right-hand side 
of Eq. (|2"3"|) which corresponds to the boundary condition 

(USD- 

In order to compute the curvature of the ionization 
front, we need to define a one-dimensional curve from 
the diffuse two-dimensional front. For this purpose, we 
use the following procedure. Let the streamer propagate 



along the x-axis, y being the transverse axis. For a given 
value of x, we locate the position of the maximum of the 
net charge density along the y-axis to get two points (due 
to mirror symmetry) of the one-dimensional curve. We 
repeat the procedure for each value x along the streamer 
length to get the complete one-dimensional curve: y(x) 
indicates the position of the maximal charge density for 
every x. The same procedure was used previously in 
Ref. (26|. We estimate the curvature n of the front by 
fitting the section of the curve y(x) with a polynomial 
x= -ny 2 /2 + 0{y 4 ). 

The enhanced field £?+( slm ) is identified with the max- 
imum of the absolute value of the electric field in the 
simulations (along the symmetry axis y = 0). 

The extraction of E~ from the simulations, E~( slm \ 
is obtained from the profile of the electric potential 
along the symmetry axis as already explained above (see 
also Fig. [2]) while its value obtained within the moving 
boundary approximation, £? - ( MBA ). is computed using 
Eq. CCH). 

The potential behind the front, <j)~ , is obtained to- 
gether with £M slm ) since the latter is given by the slope 
of the linear part of the simulated potential behind the 
front while the former is given by the intersection of the 
linear regression with the position of the tip of the front 
(the position of the discontinuity line T(t)) that here was 
chosen to be the maximum of the net charge density. 

The potential ahead of the front, tp + , is identified with 
the value of the potential at the location of the maxi- 
mum of the net charge density. We also report later in 
Fig. [5] the values of the potential at two grid points on 
our finest grid, adjacent to the location chosen to be the 
discontinuity line, r(i), of the front. 

The electron density behind the negative front, CT b s a ™\ 
is obtained from the simulations as Cback = ^(c^back) + 
c(aJend)), where Sback is the position where the net charge 
density vanishes. Such point must exist since we start 
with a neutral seed between the electrodes and thus when 
the streamer forms, positive and negative net charge den- 
sities form at the streamer edges and, consequently, the 
charge density vanishes somewhere in between. The ab- 
scissa #end is defined as Xcnd — -^max 2(x ma x 2-back); 

where x max is the position of the maximum of the 
net charge density, see Fig. [3] The quantity cr(a; enc j) 
(f (a^back)) is then the lower (upper) end of the error bars 
on the value of the electron density behind the negative 
front. This procedure gives an estimation of the interval 
of variation of a behind the front. In Fig.[3l we illustrate 
the procedure for E^ = —0.5 at time t = 490. 

The main source of errors in extracting the relevant 
quantities from the simulation data is the diffusive na- 
ture of the simulated front and hence the non- uniqueness 
in identifying the interface T(t). This uncertainty has 
no influence on the extraction of the quantities £M slm ) 

and f^ack fr' om the simulation data since those quan- 
tities are evaluated far enough behind the front where 
they are essentially constant. The error on the slope of 
the linear part of the potential behind the front, which 
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FIG. 3: (Color online) Electron and net charge density profiles 
for Eaa = —0.5 and t — 490. The position of the maximum of 
the negative net charge density, £ max , the position where the 
charge density vanishes, £back, and a; en d (see text) are also 
indicated. The distance d is equal to x max — Xhack- 



gives E ( sim ), is negligible for our purpose. The inter- 
val [x on( j, Xback] where we chose to measure cr^™' is, of 
course, somewhat arbitrary but since the electron den- 
sity is essentially constant behind the front, another pro- 
cedure would give equivalent results; only the size of the 
error bars could be slightly different. Consequently the 
errors on these two quantities are well controlled. The er- 
rors for the extracted value £'+( slm ) are also negligible for 
our purpose since this quantity is evaluated on the finest 
grid (Ax = Ay = 1/4) used in the simulations. However, 
the uncertainty about the exact position of T(t) affects 
the extraction of <p + and <$>~~ . Indeed, the location where 
we chose to evaluate cf> + and <fr~ on the potential profile 
is rather arbitrary. Moreover, the potential and the lin- 
ear regression vary significantly over short distances as 
shown in Fig. [2j Consequently, the uncertainty of the 
location of cf> + (and thus of cj>~) directly influences the 
results of the comparison between the moving boundary 
approximation and simulations. However, as explained 
in Sec. IHTCI the value [0+ - c/>~]( sim ) extracted from the 
simulation data using the procedure described above is 
an upper bound on the actual value of the potential jump. 



hanced field at the tip of the streamer. The formula (J9]) 
derived for planar fronts catches qualitatively the evolu- 
tion of the front width for a planar interface: The width 
decreases when the enhanced field increases. Moreover, 
we notice that in the present simulations (in 2D and in a 
homogeneous electric field) after some initial transients, 
the value of the enhanced electric field, up to the time of 
branching, in good approximation is given by (see Fig. 3|) 



\E+\ = 2\E n 



small corrections. 



(24) 



Consequently, the width of the front is controlled essen- 
tially by the background electric field (plus some cor- 
rections) and thus for low Eoo, where the front width 
£~ diverges, one can expect that the boundary condi- 
tions (fl"2"|) and (TiUj) and Eq. (j2"0|) will not approximate 
the simulations very well. We show below, however, that 
for Eoo > 1, our analytical approximation for the value 
of the electric field behind the front fits the simulations 
very well. 
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C. Influence of the background electric field 

We expect that the simulation results are better ap- 
proximated by the moving boundary approximation (|12[) , 
([12]) and (f2"0)) when the background electric field, E^, is 
large enough. This is so, since as mentioned above, our 
boundary conditions are derived in the regime £~n <C 1. 
The formula © and the simulations indicate that the 
width of the front is controlled by the value of the cn- 



FIG. 4: (Color online) Top: Evolution of the absolute value 
of the maximum of the electric field along the symmetry axis, 
y — 0, as a function of time for four values of the background 
electric field Eoa (B.E.F. on the figure). Middle: Evolution of 
the curvature of the tip of the front as a function of time for 
the same values of . Bottom: e as a function of time for 
the same values of E x . 

For sufficiently large E^, we expect that the relations 
(fT2")l . (fl~^ and (|20[) approximate the simulations well, and 
we also expect that the approximation improves with 
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FIG. 5: (Color online) Comparison between the simulated 
electric field behind negative ionization fronts, iJ~( slm ), and 
the values, E~^ MBA \ computed with the boundary condition 
(|12|) using curvature and enhanced field, E + , from the simu- 
lations. The comparison is performed for four values of the 
background electric field (B.E.F. in the figure). 



time. Indeed, starting from the initial neutral seed, first 
an interfacial layer forms, and then the value of the en- 
hanced field grows and approaches a plateau value given 
by (|24|) . From Eq. ([9]), this also means that the width 
of the front decreases during this time. On the other 
hand, during this process, the curvature of the front de- 
creases also (see Fig. 0]). Consequently, for a given Eoo, 
the product £~k decreases during the evolution of the 
streamer; see Fig. 01 Consequently, since our boundary 
conditions are derived for £~k -C 1, we expect better 
agreement between the moving boundary approximation 
and simulations for time and background electric fields 
large enough. 

This discussion is summarized in the lower panel of 
Fig. |H where we show that e = £~ k is a decreasing func- 
tion of time and of Eoo . 



III. RESULTS OF THE COMPARISON 



A. Testing the boundary condition for E 



Following the procedures described in Sec. IIIB1 we ex- 
tracted the values of E~( sm ^ from the simulations for 
four background electric fields: E x = —0.5,-1.0,-1.5 
and —2.0. These values are then compared with the val- 
ues, £ ,_ ( MBA ), predicted by Eq. (fT2"|) where the curvature, 
«, and the enhanced field , E + , are also obtained from 
the simulations. The results are compared in Fig. [5] The 
error bars of iJ~( slm ) are too small to be visible in the 
figure. 

The agreement between simulations and the moving 
boundary approximation is rather remarkable except for 
Eoo = —0.5. For this case, the relative errors are always 
larger than 65%, while for larger background field the er- 



rors stay always below 10 — 12%. In order to understand 
why the agreement is less good for Eoo = —0.5, we com- 
pute e from Eq. ([5]). Indeed, we recall that the moving 
boundary approximation was derived through first order 
perturbation theory in e. However, the theory, being lin- 
ear in e, does not provide an estimation for how small e 
should be. Fig.|3]shows that for Eoo = —0.5, the value of 
e stays always above 0.1. Actually from that figure, we 
can infer that for e < 0.05, the boundary condition (fT2"]) 
is accurate within 5% or less. 

However, at first sight, e seems not to be the only con- 
trol parameter. Indeed, for = —1.5 and t = 100, 
we read from Fig. 2] that e ~ 0.10 and we find that the 
relative error for E~ is about 11% (see Fig. [5]) while for 
Eoo = -0.5 and t = 490, we find that e ~ 0.11 and 
that the relative error is about 84%. This means that for 
the same value of e we get quite different relative errors 
for the values of the electric field behind the negative 
front. However for such a value of e, second order terms, 
neglected in the derivation of the moving boundary ap- 
proximation (fTO|) - (fT4")) . see Ref. [27j, could still play a 
role. For example, a coefficient associated with e 2 which 
would decrease fast enough with an increase of the en- 
hanced field may explain why second order terms are, 
in this situation, negligible for larger fields while they 
still play some role for weaker ones. Second order terms 
could also depend more significantly on the geometry of 
the streamer by involving a tangential derivative of the 
curvature. However, without deriving the second order 
theory, we cannot draw definitive conclusions on this par- 
ticular issue. 



B. Testing the relation for aback 



Using the procedure described in Sec. ITTBl we esti- 
mated the values of from the simulations for the 
same four background electric fields. These values are 
then compared with the values, c^ack ' predicted by 
Eq. (|2"0"]) . The results are compared in Fig. [51 

The simulation values and those of the moving bound- 
ary approximation agree rather well. However, the value 
of 0wv is slightly underestimated in larger fields. Nev- 
ertheless, for e < 0.05, the relative errors are about 10% 
or less. Moreover, the curvature correction improves the 
approximation of the electron density behind the front 
since the additional term is positive (see Eq. (|2"0"]) ). In 
Fig. \7\ we compare the effects of the curvature correction 
for Eoo = -1.0. 



C. Testing the boundary condition for 



Following the procedures described in Sec. IIIB1 we es- 
timated the values of [4> + — 0~]' slm - ) from the simulations 
for the same four background electric fields. These val- 
ues are then compared with the values, [cj) + — 0~]( MBA ), 
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FIG. 



(Color online) Comparison between the value of the 

and the 



simulated electron densities behind the front, c^ack' 



values, a 



(MBA) 



computed with Eq. (|20[) using curvature and 



enhanced field, E + , of the simulations. This comparison is 
performed for four values of the background electric field 
(B.E.F. in the figure). 



back 
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MBA with curvature 
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FIG. 7: (Color online) Comparison between the values of the 



simulated electron density behind the front, c 

(MBA) 



(sim) 



the val- 



'back ' 

k ' , computed with Eq. ([20)) and the values, cr^ c ^ A ^ , 
computed with k = 0. This comparison is performed for 
= -1.0. 



predicted by Eq. (fT3|) (or cquivalcntly Eq. ((23)) ). The re- 
sults are compared in Fig. [SJ The lower (upper) end of 
the error bars for the simulation results correspond to the 
value of the potential at the grid point just before (after) 
the position of the maximum of the net charge density 
on our finest grid. The size of the error bars indicates 
clearly that indeed the potential varies significantly over 
quite short distances. 

The agreement between the simulation results and the 
moving boundary approximation for the potential gap, 
Eq. (TP3"|) , is less satisfactory than the excellent agreement 
demonstrated above for the electric field and charge den- 



Sim, B.E.F. = -0.5 
MBA, B.E.F. = -0.5 
Sim, B.E.F. = -1.0 
MBA, B.E.F. = -1.0 
Sim, B.E.F. = -1.5 
MBA, B.E.F. = -1.5 
Sim, B.E.F. = -2.0 
MBA, B.E.F. 




FIG. 8: (Color online) Comparison between the jump of 
the electric potential across the interface from simulations, 
[</>+ - ci-] (sim) , and the values, [(/>+ - ^"] (MBA) , computed 
with the boundary condition (|13|) where the curvature and 
the enhanced field, E + , are also obtained from the simula- 
tions. This comparison is performed for four values of the 
background electric field (B.E.F. in the figure). 



sity, Eqs. (Q2J) and ([20)), Indeed, for e < 0.05, the relative 
error is about 20% or less while for E~ and the same 
values of e, the relative error was about 5% or less. One 
reason could simply be that the moving boundary ap- 
proximation works less well for the jump of the electric 
potential than for E~ , perhaps due to corrections asso- 
ciated with higher order terms in e. However, another 
reason is certainly that in this analysis there is one arbi- 
trariness: The precise location for evaluation of . In- 
deed, as already mentioned above, we choose the location 
of 4> + as the location of the maximum of the negative net 
charge density. Even if this is a rather natural choice, the 
actual position of cf> + , assumed by the moving boundary 
approach, could be different. However, because the po- 
tential varies significantly over short distances, see Fig. [5] 
and the size of the error bars on Fig. [8j the arbitrariness 
of the location of (f> + has certainly a direct influence on 
the comparison between the moving boundary approx- 
imation and simulations. For example, another possi- 
ble location for <fi + could be the place, x, such that the 
amount of negative charge oni<i equals the amount 
of negative charge on x > x. Since the profile of the net 
charge density is asymmetric with respect to the position 
of the maximum (see Fig. [2]), x would be located before 
the position of the maximum [x < x ma x) and the jump 
of the electric potential extracted from the data would 
be smaller since the potential and its linear regression 
are increasing functions of x. Consequently, the quantity 
[4> + — </>~]( slm ' extracted from the data using our proce- 
dure is actually an upper bound on the potential jump 
assumed in the moving boundary approach. 
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IV. CONCLUSIONS 

In this paper, we have tested the recently derived mov- 
ing boundary approximation (27j for negative ionization 
fronts on simulations of the minimal model |T])-([3]). For 
this purpose, we employed a simple geometry, charac- 
terized by mirror symmetry, and performed our analysis 
only along the symmetry axis. Future work may seek to 
perform similar analysis away from the symmetry axis. 
Our analysis confirmed the validity of two out of the three 
moving boundary conditions derived in p7j . pertaining 
to the curvature dependence of the electrostatic field and 
the charge density in the ionized region behind the propa- 
gating front. We showed that these boundary conditions 
are satisfied for slightly curved fronts, characterized by a 
small ratio between the front width £~ and the radius of 
curvature n^ 1 of the front. A third boundary condition, 
concerning the potential jump across the curved front has 
not been fully confirmed - a problem that we attribute 
to the inherent arbitrariness in extracting from simula- 
tions the appropriate potential values (corresponding to 
their value at the discontinuity line assumed by the mov- 



ing boundary approach). Further study of the range of 
validity of this condition will require the development of 
quantitative tools for such analysis. 

Finally, the usefulness of the moving boundary ap- 
proach for analytic and numerical studies of streamer 
dynamics depends crucially on its capability to describe 
front dynamics when the ratio e is not small, as could 
happen, at least in principle, along some regimes of the 
propagating front. A progress in studying this impor- 
tant question will require extension of the MBA derived 
in 27J to such regime, and comparison with numerical 



simulations along the approach developed in this paper. 
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